{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e3dd64c6-28cc-418f-ae4a-fadd4791e6b2",
   "metadata": {
    "tags": []
   },
   "source": [
    "\n",
    "# Factoring 15 with Shor's Algorithm\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbfe5758-7fa5-4596-9473-582d3197abe8",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "## Introduction\n",
    "\n",
    "The integer factorization problem [[1](#IntegerFactor)] is a famous problem in number theory: given a number $N$ which is composite, find its prime factors. The importance of the problem stems from the fact that no efficient (polynomial-time, in the number of bits needed to represent $N$) classical algorithm is known for it to this day, and much of modern day cryptography relies on this fact. In 1994, Peter Shor came up with an efficient _quantum_ algorithm for the problem [[2](#Shor94)] - providing one of the first concrete pieces of evidence for the power of quantum computers.\n",
    "\n",
    "### Shor's Algorithm\n",
    "\n",
    "Shor's algorithm consists of a classical part and a quantum subroutine. The steps of the algorithm for factoring an input number $N$, summarized from [[3](#ShorSteps)], are as follows:\n",
    "\n",
    "1. Pick a random number $1 < a < N$ that is co-prime with $N$. Co-primality can be checked by computing the GCD (greatest common divisor) of $a$ and $N$ - if it is 1 then we have found a co-prime $a$, otherwise we have found a non-trivial factor of $N$ and we are done.\n",
    "2. Find the period $r$ of the following function, using the quantum period finding algorithm (described in [[4](#PeriodFinding)]): $$f(x) = a^x \\mod N$$\n",
    "3. If $r$ is odd or $a^{r/2} = -1 \\mod N$, return to step 1 (this event can be shown to happen with probability at most $1/2$).\n",
    "4. Otherwise, $\\gcd(a^{r/2} \\pm 1, N)$ are both factors of $N$, and computing one of them yields the required result.\n",
    "\n",
    "In this demo, we will factor the number $N=15$ using Shor's algorithm, by applying the quantum subroutine (step 2) with $a=7$. This particular $a$ is chosen since it is co-prime with 15 and satisfies the conditions of step 3, providing us with a high probability of finding a factor of $N$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cc49f8e-2d6a-4a34-ae27-8e6081a10701",
   "metadata": {},
   "source": [
    "## Building the quantum period finding circuit\n",
    "\n",
    "We begin by declaring the number of qubits in the upper (counting) register the quantum subroutine uses. In our case, $N = 15$, and according to the algorithm the upper register must contain $q = \\log(Q)$ qubits for $Q$ such that $N^2 \\le Q < 2N^2$, namely $225 < Q < 450$, and therefore $q = 8$. In addition, the second register should be large enough to encode 15, hence:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8fe7f452-9fa7-484c-a824-2d521d90ee8e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:40:21.285960Z",
     "iopub.status.busy": "2024-05-07T14:40:21.285452Z",
     "iopub.status.idle": "2024-05-07T14:40:21.441875Z",
     "shell.execute_reply": "2024-05-07T14:40:21.441191Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "N = 15\n",
    "\n",
    "num_counting_qubits = int(np.ceil(np.log2(N**2)))\n",
    "num_auxilliary_qubits = int(np.ceil(np.log2(N)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21042631-a0b7-497a-9a91-2bb8e76e4562",
   "metadata": {},
   "source": [
    "We will implement a Phase Estimation [[5](#PhaseEstimation)] circuit. Each element in the circuit is a controlled operation of: $$|x\\rangle \\rightarrow |x\\cdot a^{2^i}\\mod N \\rangle $$ where $a < N$ is a number such that $\\gcd(a, N)=1$. For this demonstration we picked $a=7$. $i$ is the index of the control qubit, located in the upper register.\n",
    "\n",
    "It is quiet involved to implement these unitaries, so for this demo we will make a shortcut, and compute exactly the unitary matrix that implements the computation (which in the general case is not applicable as this pre-processing is exponential). We will do so by calculating the modular-multiplication by $a$ matrix, then using its powers.\n",
    "\n",
    "The function `unitary` is used for decomposing the unitary matrix into quantum gates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4e0570a0-efd5-4957-b75a-c7d72e5273d6",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:40:21.446692Z",
     "iopub.status.busy": "2024-05-07T14:40:21.445553Z",
     "iopub.status.idle": "2024-05-07T14:40:24.257410Z",
     "shell.execute_reply": "2024-05-07T14:40:24.256705Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from classiq import (\n",
    "    CInt,\n",
    "    Output,\n",
    "    QArray,\n",
    "    QBit,\n",
    "    X,\n",
    "    allocate,\n",
    "    control,\n",
    "    create_model,\n",
    "    hadamard_transform,\n",
    "    invert,\n",
    "    power,\n",
    "    qft,\n",
    "    qfunc,\n",
    "    repeat,\n",
    "    unitary,\n",
    ")\n",
    "\n",
    "\n",
    "def get_modular_multiplication_matrix():\n",
    "    # fmt: off\n",
    "    swap = np.array(\n",
    "        [\n",
    "            [1, 0, 0, 0],\n",
    "            [0, 0, 1, 0],\n",
    "            [0, 1, 0, 0],\n",
    "            [0, 0, 0, 1]\n",
    "        ],\n",
    "        dtype=complex\n",
    "    )\n",
    "    # fmt: on\n",
    "    swap32 = np.kron(np.identity(4), swap)\n",
    "    swap21 = np.kron(np.kron(np.identity(2), swap), np.identity(2))\n",
    "    swap10 = np.kron(swap, np.identity(4))\n",
    "    x = np.array([[0, 1], [1, 0]])\n",
    "    x_all = np.kron(np.kron(x, x), np.kron(x, x))\n",
    "    u = x_all @ swap10 @ swap21 @ swap32\n",
    "    return u\n",
    "\n",
    "\n",
    "MODULAR_MUL_UNITARY = get_modular_multiplication_matrix().real.tolist()\n",
    "\n",
    "\n",
    "@qfunc\n",
    "def modular_exponentiation(\n",
    "    exponent: CInt, target: QArray[QBit, num_auxilliary_qubits]\n",
    ") -> None:\n",
    "    power(2**exponent, lambda: unitary(elements=MODULAR_MUL_UNITARY, target=target))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "545bb4fd-ff8f-4ff8-b86d-d89808d91abb",
   "metadata": {},
   "source": [
    "### Building the complete circuit"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a97036e1-87cf-4506-ae61-dd15588ec247",
   "metadata": {},
   "source": [
    "At the first layer of the quantum circuit, we prepare the equal superposition state in the top (counting) register, and prepare the $|1\\rangle$ state in the bottom (auxiliary) register."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05072f38-8ffd-482c-9af0-8aa4b404c2f1",
   "metadata": {},
   "source": [
    "We then apply the second layer of the circuit, which consists of the controlled $U^{2^i}$ gates. \n",
    "Lastly, we apply an inverse QFT on the counting register, to get the period."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b2c2aeac-2049-45fe-a529-f98aed009b37",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:40:24.260542Z",
     "iopub.status.busy": "2024-05-07T14:40:24.259953Z",
     "iopub.status.idle": "2024-05-07T14:40:24.264624Z",
     "shell.execute_reply": "2024-05-07T14:40:24.263949Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "@qfunc\n",
    "def period_finding(\n",
    "    qv_counting: Output[QArray[QBit, num_counting_qubits]],\n",
    "    qv_auxilliary: Output[QArray[QBit, num_auxilliary_qubits]],\n",
    ") -> None:\n",
    "    # start with a hadamard transform in the counting register\n",
    "    allocate(num_counting_qubits, qv_counting)\n",
    "    hadamard_transform(qv_counting)\n",
    "\n",
    "    # Prepare the |1> state on the lower register\n",
    "    allocate(num_auxilliary_qubits, qv_auxilliary)\n",
    "    X(qv_auxilliary[0])\n",
    "\n",
    "    # Apply the contolled modular-exponentiations using each of the counting qubits\n",
    "    repeat(\n",
    "        count=num_auxilliary_qubits,\n",
    "        iteration=lambda index: control(\n",
    "            ctrl=qv_counting[index],\n",
    "            operand=lambda: modular_exponentiation(index, qv_auxilliary),\n",
    "        ),\n",
    "    )  # ! not working with qv[a:]\n",
    "\n",
    "    # Lastly, apply an inverse QFT\n",
    "    invert(lambda: qft(qv_counting))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "579c9843-907e-4454-a92c-6a0a04d0615c",
   "metadata": {},
   "source": [
    "### Quantum entry point\n",
    "In order to synthesize the circuit, we define a quantum `main` function. As are we only interested in the output of the counting register, we only define it in the signature of the function.\n",
    "\n",
    "Next, we translate it to qmod using the `create_model`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4d7da11c-adfc-4ac8-9277-0c97fd65539e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:40:24.266926Z",
     "iopub.status.busy": "2024-05-07T14:40:24.266735Z",
     "iopub.status.idle": "2024-05-07T14:40:24.285019Z",
     "shell.execute_reply": "2024-05-07T14:40:24.284422Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "@qfunc\n",
    "def main(qv_counting: Output[QArray[QBit, num_counting_qubits]]) -> None:\n",
    "    qv_auxilliary = QArray(\"qv_auxilliary\")\n",
    "    period_finding(qv_counting, qv_auxilliary)\n",
    "\n",
    "\n",
    "qmod = create_model(main)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "3adf6690-8f4c-495e-9ddc-f4702f1efcdd",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:40:24.287286Z",
     "iopub.status.busy": "2024-05-07T14:40:24.287108Z",
     "iopub.status.idle": "2024-05-07T14:40:24.307133Z",
     "shell.execute_reply": "2024-05-07T14:40:24.306559Z"
    }
   },
   "outputs": [],
   "source": [
    "from classiq import write_qmod\n",
    "\n",
    "write_qmod(qmod, \"shor\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08e1a6a0-137d-4c49-a215-97daa2197f5c",
   "metadata": {},
   "source": [
    "We now send the model to the synthesis engine, taking a few seconds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d1eb4133-896d-4e2d-8253-83eb68bf467d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:40:24.310100Z",
     "iopub.status.busy": "2024-05-07T14:40:24.309720Z",
     "iopub.status.idle": "2024-05-07T14:40:37.844547Z",
     "shell.execute_reply": "2024-05-07T14:40:37.843782Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from classiq import synthesize\n",
    "\n",
    "qprog = synthesize(qmod)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47ce0b6d-839a-4b4c-82cd-daf4b9b5be14",
   "metadata": {},
   "source": [
    "We can now view the circuit and its depth:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "08ce2565-5197-4c7a-b7c5-1dad94e4e46a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:40:37.850425Z",
     "iopub.status.busy": "2024-05-07T14:40:37.850155Z",
     "iopub.status.idle": "2024-05-07T14:40:38.022304Z",
     "shell.execute_reply": "2024-05-07T14:40:38.021691Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Opening: https://platform.classiq.io/circuit/86e59cb0-aa57-4c3f-933a-c33a0e05bbdc?version=0.41.0.dev39%2B79c8fd0855\n"
     ]
    }
   ],
   "source": [
    "from classiq import show\n",
    "\n",
    "show(qprog)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fc7fe66-c1ca-46b0-9272-cfcc3c6e73f3",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Executing the circuit\n",
    "\n",
    "Now, we turn to executing the circuit above, using the simulator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "3f2ba37b-f01b-48f2-85a0-ccae575d3de5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:40:38.026756Z",
     "iopub.status.busy": "2024-05-07T14:40:38.025681Z",
     "iopub.status.idle": "2024-05-07T14:40:39.417985Z",
     "shell.execute_reply": "2024-05-07T14:40:39.417368Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from classiq import execute\n",
    "\n",
    "results = execute(qprog).result()\n",
    "res = results[0].value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "a7c5710d-f360-4500-b489-0986b7178938",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:40:39.420639Z",
     "iopub.status.busy": "2024-05-07T14:40:39.420272Z",
     "iopub.status.idle": "2024-05-07T14:40:39.423609Z",
     "shell.execute_reply": "2024-05-07T14:40:39.423037Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import collections\n",
    "\n",
    "hist_counting_qubits = collections.defaultdict(int)\n",
    "for key, value in res.counts.items():\n",
    "    hist_counting_qubits[key] += value"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c351af0c-696b-4ba9-b671-f9759c1ca387",
   "metadata": {},
   "source": [
    "Plotting the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "b0f701d4-69ff-4263-aba4-94a234be8f4e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T14:40:39.426056Z",
     "iopub.status.busy": "2024-05-07T14:40:39.425714Z",
     "iopub.status.idle": "2024-05-07T14:40:39.553081Z",
     "shell.execute_reply": "2024-05-07T14:40:39.552441Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<BarContainer object of 4 artists>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAfdklEQVR4nO3df2xV9f3H8Vd/0Esr3NaCvZfOgvgTEAQGWu5EwqSh/AjqbDJhTNEQmK6YaB2ybgiCy7eOGTWSCnNRmImIkihGZJ1YlOosKF2Rn2uE4cDBLRNsLzApLf18/1g48UqhFPrj3evzkZyEe87nnvu59+SUZ+69p41zzjkBAAAYEt/REwAAAPguAgUAAJhDoAAAAHMIFAAAYA6BAgAAzCFQAACAOQQKAAAwh0ABAADmJHb0BC5EY2OjDhw4oO7duysuLq6jpwMAAM6Dc05Hjx5VZmam4uPP/R5JpwyUAwcOKCsrq6OnAQAALsD+/ft1+eWXn3NMpwyU7t27S/rfE/T7/R08GwAAcD4ikYiysrK8/8fPpVMGyumPdfx+P4ECAEAncz5fz+BLsgAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5iR09AQC4WFf8+p2OnsL31hdPTuzoKSBG8Q4KAAAwh0ABAADmECgAAMAcAgUAAJhDoAAAAHMIFAAAYA6BAgAAzCFQAACAOQQKAAAwh0ABAADmECgAAMAcAgUAAJhDoAAAAHMIFAAAYA6BAgAAzCFQAACAOQQKAAAwh0ABAADmECgAAMAcAgUAAJhDoAAAAHMSO3oCAACczRW/fqejp/C99cWTEzv08XkHBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAc1oUKEVFRbrxxhvVvXt3ZWRk6I477lBVVVXUmNGjRysuLi5quf/++6PG7Nu3TxMnTlRKSooyMjI0e/ZsNTQ0XPyzAQAAMaFFf4tnw4YNys/P14033qiGhgb95je/0dixY7Vz505dcskl3rgZM2Zo4cKF3u2UlBTv36dOndLEiRMVDAb18ccf6+DBg7rnnnvUpUsX/d///V8rPCUAANDZtShQSkpKom4vX75cGRkZqqio0KhRo7z1KSkpCgaDTe7j3Xff1c6dO/Xee+8pEAhoyJAheuKJJzRnzhw9/vjjSkpKuoCnAQAAYslFfQeltrZWkpSenh61/pVXXlHPnj01cOBAFRYW6r///a+3rby8XIMGDVIgEPDW5ebmKhKJaMeOHU0+Tl1dnSKRSNQCAABiV4veQfm2xsZGPfTQQ7r55ps1cOBAb/3PfvYz9enTR5mZmdq6davmzJmjqqoqvfHGG5KkcDgcFSeSvNvhcLjJxyoqKtKCBQsudKoAAKCTueBAyc/P1/bt2/XRRx9FrZ85c6b370GDBqlXr14aM2aM9uzZo6uuuuqCHquwsFAFBQXe7UgkoqysrAubOAAAMO+CPuKZNWuW1qxZo/fff1+XX375OcdmZ2dLknbv3i1JCgaDqq6ujhpz+vbZvrfi8/nk9/ujFgAAELtaFCjOOc2aNUtvvvmm1q9fr759+zZ7ny1btkiSevXqJUkKhULatm2bDh065I1Zt26d/H6/BgwY0JLpAACAGNWij3jy8/O1YsUKvfXWW+revbv3nZHU1FQlJydrz549WrFihSZMmKAePXpo69atevjhhzVq1CjdcMMNkqSxY8dqwIABuvvuu7Vo0SKFw2HNnTtX+fn58vl8rf8MAQBAp9Oid1CWLFmi2tpajR49Wr169fKW1157TZKUlJSk9957T2PHjlW/fv30yCOPKC8vT2+//ba3j4SEBK1Zs0YJCQkKhUL6+c9/rnvuuSfq96YAAIDvtxa9g+KcO+f2rKwsbdiwodn99OnTR2vXrm3JQwMAgO8R/hYPAAAwh0ABAADmECgAAMAcAgUAAJhDoAAAAHMIFAAAYA6BAgAAzCFQAACAOQQKAAAwp0W/SRbozK749TsdPYXvrS+enNjRUwDQyfAOCgAAMIdAAQAA5hAoAADAHAIFAACYQ6AAAABzCBQAAGAOgQIAAMwhUAAAgDkECgAAMIdAAQAA5hAoAADAHAIFAACYQ6AAAABzCBQAAGAOgQIAAMwhUAAAgDkECgAAMIdAAQAA5iR29AQsuuLX73T0FL63vnhyYkdPAQBgAO+gAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwJwWBUpRUZFuvPFGde/eXRkZGbrjjjtUVVUVNebEiRPKz89Xjx491K1bN+Xl5am6ujpqzL59+zRx4kSlpKQoIyNDs2fPVkNDw8U/GwAAEBNaFCgbNmxQfn6+Nm7cqHXr1qm+vl5jx47V8ePHvTEPP/yw3n77ba1atUobNmzQgQMHdOedd3rbT506pYkTJ+rkyZP6+OOP9ec//1nLly/XvHnzWu9ZAQCATi2xJYNLSkqibi9fvlwZGRmqqKjQqFGjVFtbqxdffFErVqzQrbfeKklatmyZ+vfvr40bN2rEiBF69913tXPnTr333nsKBAIaMmSInnjiCc2ZM0ePP/64kpKSWu/ZAQCATumivoNSW1srSUpPT5ckVVRUqL6+Xjk5Od6Yfv36qXfv3iovL5cklZeXa9CgQQoEAt6Y3NxcRSIR7dixo8nHqaurUyQSiVoAAEDsuuBAaWxs1EMPPaSbb75ZAwcOlCSFw2ElJSUpLS0tamwgEFA4HPbGfDtOTm8/va0pRUVFSk1N9ZasrKwLnTYAAOgELjhQ8vPztX37dq1cubI159OkwsJC1dbWesv+/fvb/DEBAEDHadF3UE6bNWuW1qxZo7KyMl1++eXe+mAwqJMnT6qmpibqXZTq6moFg0FvzCeffBK1v9NX+Zwe810+n08+n+9CpgoAADqhFr2D4pzTrFmz9Oabb2r9+vXq27dv1PZhw4apS5cuKi0t9dZVVVVp3759CoVCkqRQKKRt27bp0KFD3ph169bJ7/drwIABF/NcAABAjGjROyj5+flasWKF3nrrLXXv3t37zkhqaqqSk5OVmpqq6dOnq6CgQOnp6fL7/XrwwQcVCoU0YsQISdLYsWM1YMAA3X333Vq0aJHC4bDmzp2r/Px83iUBAACSWhgoS5YskSSNHj06av2yZct07733SpKeeeYZxcfHKy8vT3V1dcrNzdXzzz/vjU1ISNCaNWv0wAMPKBQK6ZJLLtG0adO0cOHCi3smAAAgZrQoUJxzzY7p2rWriouLVVxcfNYxffr00dq1a1vy0AAA4HuEv8UDAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOa0OFDKyso0adIkZWZmKi4uTqtXr47afu+99youLi5qGTduXNSYI0eOaOrUqfL7/UpLS9P06dN17Nixi3oiAAAgdrQ4UI4fP67BgweruLj4rGPGjRungwcPesurr74atX3q1KnasWOH1q1bpzVr1qisrEwzZ85s+ewBAEBMSmzpHcaPH6/x48efc4zP51MwGGxy265du1RSUqJPP/1Uw4cPlyQtXrxYEyZM0FNPPaXMzMyWTgkAAMSYNvkOygcffKCMjAxdd911euCBB3T48GFvW3l5udLS0rw4kaScnBzFx8dr06ZNTe6vrq5OkUgkagEAALGr1QNl3Lhxevnll1VaWqrf//732rBhg8aPH69Tp05JksLhsDIyMqLuk5iYqPT0dIXD4Sb3WVRUpNTUVG/Jyspq7WkDAABDWvwRT3MmT57s/XvQoEG64YYbdNVVV+mDDz7QmDFjLmifhYWFKigo8G5HIhEiBQCAGNbmlxlfeeWV6tmzp3bv3i1JCgaDOnToUNSYhoYGHTly5KzfW/H5fPL7/VELAACIXW0eKF9++aUOHz6sXr16SZJCoZBqampUUVHhjVm/fr0aGxuVnZ3d1tMBAACdQIs/4jl27Jj3bogk7d27V1u2bFF6errS09O1YMEC5eXlKRgMas+ePXr00Ud19dVXKzc3V5LUv39/jRs3TjNmzNDSpUtVX1+vWbNmafLkyVzBAwAAJF3AOyibN2/W0KFDNXToUElSQUGBhg4dqnnz5ikhIUFbt27VbbfdpmuvvVbTp0/XsGHD9OGHH8rn83n7eOWVV9SvXz+NGTNGEyZM0MiRI/XCCy+03rMCAACdWovfQRk9erScc2fd/te//rXZfaSnp2vFihUtfWgAAPA9wd/iAQAA5hAoAADAHAIFAACYQ6AAAABzCBQAAGAOgQIAAMwhUAAAgDkECgAAMIdAAQAA5hAoAADAHAIFAACYQ6AAAABzCBQAAGAOgQIAAMwhUAAAgDkECgAAMIdAAQAA5hAoAADAHAIFAACYQ6AAAABzCBQAAGAOgQIAAMwhUAAAgDkECgAAMIdAAQAA5hAoAADAHAIFAACYQ6AAAABzCBQAAGAOgQIAAMwhUAAAgDkECgAAMIdAAQAA5hAoAADAHAIFAACYQ6AAAABzCBQAAGAOgQIAAMwhUAAAgDkECgAAMIdAAQAA5hAoAADAHAIFAACYQ6AAAABzCBQAAGAOgQIAAMwhUAAAgDkECgAAMIdAAQAA5hAoAADAHAIFAACYQ6AAAABzWhwoZWVlmjRpkjIzMxUXF6fVq1dHbXfOad68eerVq5eSk5OVk5Ojzz//PGrMkSNHNHXqVPn9fqWlpWn69Ok6duzYRT0RAAAQO1ocKMePH9fgwYNVXFzc5PZFixbpueee09KlS7Vp0yZdcsklys3N1YkTJ7wxU6dO1Y4dO7Ru3TqtWbNGZWVlmjlz5oU/CwAAEFMSW3qH8ePHa/z48U1uc87p2Wef1dy5c3X77bdLkl5++WUFAgGtXr1akydP1q5du1RSUqJPP/1Uw4cPlyQtXrxYEyZM0FNPPaXMzMyLeDoAACAWtOp3UPbu3atwOKycnBxvXWpqqrKzs1VeXi5JKi8vV1pamhcnkpSTk6P4+Hht2rSpyf3W1dUpEolELQAAIHa1aqCEw2FJUiAQiFofCAS8beFwWBkZGVHbExMTlZ6e7o35rqKiIqWmpnpLVlZWa04bAAAY0ymu4iksLFRtba237N+/v6OnBAAA2lCrBkowGJQkVVdXR62vrq72tgWDQR06dChqe0NDg44cOeKN+S6fzye/3x+1AACA2NWqgdK3b18Fg0GVlpZ66yKRiDZt2qRQKCRJCoVCqqmpUUVFhTdm/fr1amxsVHZ2dmtOBwAAdFItvorn2LFj2r17t3d779692rJli9LT09W7d2899NBD+t3vfqdrrrlGffv21WOPPabMzEzdcccdkqT+/ftr3LhxmjFjhpYuXar6+nrNmjVLkydP5goeAAAg6QICZfPmzfrxj3/s3S4oKJAkTZs2TcuXL9ejjz6q48ePa+bMmaqpqdHIkSNVUlKirl27evd55ZVXNGvWLI0ZM0bx8fHKy8vTc8891wpPBwAAxIIWB8ro0aPlnDvr9ri4OC1cuFALFy4865j09HStWLGipQ8NAAC+JzrFVTwAAOD7hUABAADmECgAAMAcAgUAAJhDoAAAAHMIFAAAYA6BAgAAzCFQAACAOQQKAAAwh0ABAADmECgAAMAcAgUAAJhDoAAAAHMIFAAAYA6BAgAAzCFQAACAOQQKAAAwh0ABAADmECgAAMAcAgUAAJhDoAAAAHMIFAAAYA6BAgAAzCFQAACAOQQKAAAwh0ABAADmECgAAMAcAgUAAJhDoAAAAHMIFAAAYA6BAgAAzCFQAACAOQQKAAAwh0ABAADmECgAAMAcAgUAAJhDoAAAAHMIFAAAYA6BAgAAzCFQAACAOQQKAAAwh0ABAADmECgAAMAcAgUAAJhDoAAAAHMIFAAAYA6BAgAAzCFQAACAOQQKAAAwh0ABAADmECgAAMAcAgUAAJjT6oHy+OOPKy4uLmrp16+ft/3EiRPKz89Xjx491K1bN+Xl5am6urq1pwEAADqxNnkH5frrr9fBgwe95aOPPvK2Pfzww3r77be1atUqbdiwQQcOHNCdd97ZFtMAAACdVGKb7DQxUcFg8Iz1tbW1evHFF7VixQrdeuutkqRly5apf//+2rhxo0aMGNEW0wEAAJ1Mm7yD8vnnnyszM1NXXnmlpk6dqn379kmSKioqVF9fr5ycHG9sv3791Lt3b5WXl591f3V1dYpEIlELAACIXa0eKNnZ2Vq+fLlKSkq0ZMkS7d27V7fccouOHj2qcDispKQkpaWlRd0nEAgoHA6fdZ9FRUVKTU31lqysrNaeNgAAMKTVP+IZP3689+8bbrhB2dnZ6tOnj15//XUlJydf0D4LCwtVUFDg3Y5EIkQKAAAxrM0vM05LS9O1116r3bt3KxgM6uTJk6qpqYkaU11d3eR3Vk7z+Xzy+/1RCwAAiF1tHijHjh3Tnj171KtXLw0bNkxdunRRaWmpt72qqkr79u1TKBRq66kAAIBOotU/4vnVr36lSZMmqU+fPjpw4IDmz5+vhIQETZkyRampqZo+fboKCgqUnp4uv9+vBx98UKFQiCt4AACAp9UD5csvv9SUKVN0+PBhXXbZZRo5cqQ2btyoyy67TJL0zDPPKD4+Xnl5eaqrq1Nubq6ef/751p4GAADoxFo9UFauXHnO7V27dlVxcbGKi4tb+6EBAECM4G/xAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5HRooxcXFuuKKK9S1a1dlZ2frk08+6cjpAAAAIzosUF577TUVFBRo/vz5+vvf/67BgwcrNzdXhw4d6qgpAQAAIzosUJ5++mnNmDFD9913nwYMGKClS5cqJSVFL730UkdNCQAAGJHYEQ968uRJVVRUqLCw0FsXHx+vnJwclZeXnzG+rq5OdXV13u3a2lpJUiQSaZP5Ndb9t032i+a11TGVOK4dqS2Pq8Sx7Ugc29jVFsf29D6dc82O7ZBA+eqrr3Tq1CkFAoGo9YFAQP/4xz/OGF9UVKQFCxacsT4rK6vN5oiOkfpsR88AbYHjGrs4trGrLY/t0aNHlZqaes4xHRIoLVVYWKiCggLvdmNjo44cOaIePXooLi6uA2dmSyQSUVZWlvbv3y+/39/R00Er4tjGJo5r7OLYNs05p6NHjyozM7PZsR0SKD179lRCQoKqq6uj1ldXVysYDJ4x3ufzyefzRa1LS0tryyl2an6/nxMiRnFsYxPHNXZxbM/U3Dsnp3XIl2STkpI0bNgwlZaWeusaGxtVWlqqUCjUEVMCAACGdNhHPAUFBZo2bZqGDx+um266Sc8++6yOHz+u++67r6OmBAAAjOiwQLnrrrv0n//8R/PmzVM4HNaQIUNUUlJyxhdncf58Pp/mz59/xsdh6Pw4trGJ4xq7OLYXL86dz7U+AAAA7Yi/xQMAAMwhUAAAgDkECgAAMIdAAQAA5hAoraisrEyTJk1SZmam4uLitHr16qjtzjnNmzdPvXr1UnJysnJycvT5559HjTly5IimTp0qv9+vtLQ0TZ8+XceOHYsas3XrVt1yyy3q2rWrsrKytGjRojPmsmrVKvXr109du3bVoEGDtHbt2g6bS6woLi7WFVdcoa5duyo7O1uffPKJt+3EiRPKz89Xjx491K1bN+Xl5Z3xiwj37duniRMnKiUlRRkZGZo9e7YaGhqixnzwwQf64Q9/KJ/Pp6uvvlrLly9v0Tzaey6x4Fyv5wsvvKDRo0fL7/crLi5ONTU1Z9yfc9am5n4ev/HGGxo7dqz3G8m3bNlyxj5i8bzuVBxazdq1a91vf/tb98YbbzhJ7s0334za/uSTT7rU1FS3evVq99lnn7nbbrvN9e3b133zzTfemHHjxrnBgwe7jRs3ug8//NBdffXVbsqUKd722tpaFwgE3NSpU9327dvdq6++6pKTk90f//hHb8zf/vY3l5CQ4BYtWuR27tzp5s6d67p06eK2bdvW7nOJFStXrnRJSUnupZdecjt27HAzZsxwaWlprrq62jnn3P333++ysrJcaWmp27x5sxsxYoT70Y9+5N2/oaHBDRw40OXk5LjKykq3du1a17NnT1dYWOiN+ec//+lSUlJcQUGB27lzp1u8eLFLSEhwJSUl5z2P9pxLLGju9XzmmWdcUVGRKyoqcpLc119/fcY+OGdtau7n8csvv+wWLFjg/vSnPzlJrrKy8ox9xNp53dkQKG3kuydEY2OjCwaD7g9/+IO3rqamxvl8Pvfqq68655zbuXOnk+Q+/fRTb8xf/vIXFxcX5/79738755x7/vnn3aWXXurq6uq8MXPmzHHXXXedd/unP/2pmzhxYtR8srOz3S9+8Yt2n0usuOmmm1x+fr53+9SpUy4zM9MVFRW5mpoa16VLF7dq1Spv+65du5wkV15e7pz73w/L+Ph4Fw6HvTFLlixxfr/fe/0effRRd/3110c97l133eVyc3PPax7OuXadSyxo7vU87f33328yUDhnO4emAuW0vXv3NhkosXhedzZ8xNNO9u7dq3A4rJycHG9damqqsrOzVV5eLkkqLy9XWlqahg8f7o3JyclRfHy8Nm3a5I0ZNWqUkpKSvDG5ubmqqqrS119/7Y359uOcHnP6cdpzLrHg5MmTqqioiHq94uPjlZOTo/LyclVUVKi+vj5qe79+/dS7d++o13PQoEFRv4gwNzdXkUhEO3bs8Mac67g1Nw9J7TaXWHA+r2dzOGdjVyye150NgdJOwuGwJJ3xm3IDgYC3LRwOKyMjI2p7YmKi0tPTo8Y0tY9vP8bZxnx7e3vNJRZ89dVXOnXq1Flfr3A4rKSkpDP+gOV3X88LPW6RSETffPNNs/M4vY/2mEssOJ/Xszmcs7ErFs/rzoZAAQAA5hAo7SQYDErSGd+6rq6u9rYFg0EdOnQoantDQ4OOHDkSNaapfXz7Mc425tvb22susaBnz55KSEg46+sVDAZ18uTJM67w+O7reaHHze/3Kzk5udl5nN5He8wlFpzP69kcztnYFYvndWdDoLSTvn37KhgMqrS01FsXiUS0adMmhUIhSVIoFFJNTY0qKiq8MevXr1djY6Oys7O9MWVlZaqvr/fGrFu3Ttddd50uvfRSb8y3H+f0mNOP055ziQVJSUkaNmxY1OvV2Nio0tJShUIhDRs2TF26dInaXlVVpX379kW9ntu2bYv6D2TdunXy+/0aMGCAN+Zcx625eUhqt7nEgvN5PZvDORu7YvG87nQ6+lu6seTo0aOusrLSVVZWOknu6aefdpWVle5f//qXc+5/lwmmpaW5t956y23dutXdfvvtTV4mOHToULdp0yb30UcfuWuuuSbqMsGamhoXCATc3Xff7bZv3+5WrlzpUlJSzrhkMTEx0T311FNu165dbv78+U1estgec4kVK1eudD6fzy1fvtzt3LnTzZw506WlpXnfmL///vtd79693fr1693mzZtdKBRyoVDIu//pSwDHjh3rtmzZ4kpKStxll13W5OWIs2fPdrt27XLFxcVNXo54rnm051xiQXOv58GDB11lZaV3KWpZWZmrrKx0hw8f9vbBOWtTcz+PDx8+7CorK90777zjJLmVK1e6yspKd/DgQW8fsXZedzYESis6fSnid5dp06Y55/53qeBjjz3mAoGA8/l8bsyYMa6qqipqH4cPH3ZTpkxx3bp1c36/3913333u6NGjUWM+++wzN3LkSOfz+dwPfvAD9+STT54xl9dff91de+21LikpyV1//fXunXfeidrennOJFYsXL3a9e/d2SUlJ7qabbnIbN270tn3zzTful7/8pbv00ktdSkqK+8lPfhL1g84557744gs3fvx4l5yc7Hr27OkeeeQRV19fHzXm/fffd0OGDHFJSUnuyiuvdMuWLWvRPNp7LrHgXK/n/Pnzmzynv/1acM7a1NzP42XLljW5ff78+d4+YvG87kzinHOu/d6vAQAAaB7fQQEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAc/4f2ln94fAp6fUAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.bar(hist_counting_qubits.keys(), hist_counting_qubits.values())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "993ec133-5185-4aec-a396-b0cb6762e9bb",
   "metadata": {},
   "source": [
    "We obtained 4 results $y$ from the circuit, each with probability roughly $1/4$: $0, 64, 128$ and $192$. Dividing by $Q = 256$ we obtain 4 reduced fractions: $0, 1/4, 1/2$ and $3/4$, with two of them having the correct period $r=4$ in the denominator. With this period, we can compute the factors of $N = 15$: $\\gcd(a^{r/2} \\pm 1, N) = \\gcd(7^2 \\pm 1, 15) = 3, 5$.\n",
    "\n",
    "## References\n",
    "\n",
    "<a id='IntegerFactor'>[1]</a>: [Integer Factorization (Wikipedia)](https://en.wikipedia.org/wiki/Integer_factorization)\n",
    "\n",
    "<a id='Shor94'>[2]</a>: [Shor, Peter W. \"Algorithms for quantum computation: discrete logarithms and factoring.\" Proceedings 35th annual symposium on foundations of computer science. Ieee, 1994.](https://ieeexplore.ieee.org/abstract/document/365700)\n",
    "\n",
    "<a id='ShorSteps'>[3]</a>: [Shor's Algorithm Procedure (Wikipedia)](https://en.wikipedia.org/wiki/Shor%27s_algorithm#Procedure)\n",
    "\n",
    "<a id='PeriodFinding'>[4]</a>: [Quantum Period Finding (Wikipedia)](https://en.wikipedia.org/wiki/Shor%27s_algorithm#Quantum_part:_period-finding_subroutine)\n",
    "\n",
    "<a id='PhaseEstimation'>[5]</a>: [Quantum Phase Estimation (Wikipedia)](https://en.wikipedia.org/wiki/Quantum_phase_estimation_algorithm)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
